Enhancement of Wigner crystallization in quasi low-dimensional solids. 
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The crystallization of electrons in quasi low-dimensional solids is studied in a model which retains 
the full three-dimensional nature of the Coulomb interactions. We show that restricting the electron 
motion to layers (or chains) gives rise to a rich sequence of structural transitions upon varying the 
particle density. In addition, the concurrence of low-dimensional electron motion and isotropic 
Coulomb interactions leads to a sizeable stabilization of the Wigner crystal, which could be one of 
the mechanisms at the origin of the charge ordered phases frequently observed in such compounds. 



I. INTRODUCTION 

Despite being a well established concept in the 
physics of interacting electrons, direct evidence of Wigner 
crystallization^ has been reported unambiguously in only 
a limited number of systems, namely electrons at the 
surface of liquid helium; 2 ^** and in semiconductor het- 
erostructures of extreme purity^ In both cases, a two 
dimensional electron gas (2DEG) is realized at an inter- 
face between two media, for which the jellium model of 
the homogeneous electron gas constitutes a good approx- 
imation. 

Alternatively, the charge ordering phenomena observed 
at low temperatures in a number of solids have been ei- 
ther interpreted as some form of Wigner crystallization, 
or ascribed to the presence of long-ranged Coulomb in- 
teractions. These include the one-dimensional organic 
salts TTF-TCNQ, 6 (TMTTF) 2 X, 7 - 8 (DI-DCNQI) 2 Ag^ 
the ladder cuprate compounds Sri4Cu2404i^Siii and 
chain compounds Nai+zCuO^r^ as wei l as the lay- 
ered superconducting cuprate a 14 i 15 and possibly the two- 
dimensional BEDT-TTF organic salts^ For such sys- 
tems, the jellium model is a priori a rather crude mod- 
elization, and the concept of Wigner crystallization must 
be generalized to account for other competing effects such 
as the periodic potential of the underlying lattice, chem- 
ical impurities, structural defects, magnetic interactions, 
etc. In narrow band solids, for instance, the interplay 
with the host lattice of ions can strongly affect the charge 
ordering pattern especially at highly commensurate band 
fillingspiii& Nevertheless, when the radius of localization 
of the particles is larger than the typical ion-ion distance, 
the host lattice can be replaced to a good accuracy by 
an effective continuous medium, restoring de facto the 
validity of the jellium model>iii& 

Setting aside the important problem of the commen- 
surability with the host lattice, and neglecting disorder 
and other effects that can certainly play a role in the 
compounds under study, we come to the following ob- 
servation: a common feature shared by the experimen- 
tal systems listed above is that they are all quasi low- 
dimensional solids, i.e. they are bulk three-dimensional 
(3D) compounds where the transfer integrals between dif- 



ferent chemical units are so anisotropic that the carrier 
motion is effectively restricted to two-dimensional (2D) 
atomic layers, or one-dimensional (ID) chains. Yet, the 
Coulomb forces retain their three-dimensional character, 
being long-ranged and isotropic. In such systems, inter- 
layer (inter-chain) interactions cannot be neglected, lead- 
ing eventually to a full three-dimensional ordering of the 
chargesiiiii^ This suggests why quasi low-dimensional 
solids are a particularly favorable ground for the obser- 
vation of Wigner crystallization: the electron-electron 
interactions have the same behavior as in bulk three- 
dimensional systems, but the kinetic part is strongly re- 
duced by the effective lowering of dimensionality. Re- 
minding that a Wigner crystal arises from the competi- 
tion between potential and kinetic energy, this results in 
a sizeable stabilization of the crystal as compared with 
the usual 3D case^ 

A similar conclusion is reached by observing that, even 
compared to purely low-dimensional systems such as the 
2DEG mentioned above, the Wigner crystal phase could 
be stabilized in quasi low-dimensional solids due to the 
presence of additional interlayer interactions. This topic 
has been analyzed in the literature in the framework of 
bilayer quantum wells, i.e. constituted of two coupled 
2D electron systems, where it has been shown that, de- 
pending on the strength of the interlayer forces, the or- 
dering pattern can differ from the hexagonal structure 
expected in a single layer .isLSl More importantly, it was 
founo— fi±^ that at interlayer separations comparable 
with the mean interparticle distance, the melting density 
is raised by a factor of 3 with respect to the pure 2D case, 
which makes a factor as large as 10 2 when appropriately 
scaled to the 3D situation. 

In this work, we model quasi two-dimensional (one- 
dimensional) systems as periodic arrays of conducting 
layers (wires) embedded in a three-dimensional bulk ma- 
terial, where the electrons interact through isotropic 
long-range Coulomb forces. We show that, upon vary- 
ing the particle density or the interlayer (interwire) sep- 
aration, the Wigner crystal undergoes several structural 
transitions in order to minimize its energy compatibly 
with the given geometrical constraints. We then give a 
semi-quantitative estimate of the melting density for the 
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different structures previously identified, based on the 
Lindemann criterion, which confirms the stabilization of 
the crystallized phase expected from general grounds. 

The paper is organized as follows: In Section II, we 
introduce a model for the crystallization of electrons in 
an anisotropic environment and the method for calculat- 
ing the crystal energy in the harmonic approximation, 
which includes the classical Madelung energy and the 
zero-point vibrational energy of the collective excitations. 
This is applied to the case of quasi two-dimensional sys- 
tems, for which the structural/melting phase diagram is 
determined. The validity of the present approximation 
scheme is checked at the end of Section II by analyz- 
ing a system of two coupled layers, for which our results 
compare positively with the numerical results available 
in the literature. An analogous discussion for quasi one- 
dimensional systems is reported in Section III, by treat- 
ing explicitely the case where the conducting chains form 
a square array. The main results are summarized in sec- 
tion IV. 



II. WIGNER CRYSTALLIZATION IN LAYERED 
SOLIDS 

A. Model and approximations 

Let us consider a system of electrons (or holes) of den- 
sity n = (47rr^/3) _1 in a strongly anisotropic environ- 
ment, such that the particle motion is constrained to 
equally spaced atomic layers (at distance d), but remains 
isotropic within the layers. To ensure charge neutral- 
ity, we assume a uniform 3D compensating background 
of opposite charge. The hamiltonian for N crystallized 
particles in a volume V is given by: 



H = NE M + 
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is the Madelung energy of the given lattice structure (in 
the thermodynamic limit, N, V — > oo, boundary effects 
are negligible and all particles become equivalent). The 
second term is the two-dimensional kinetic energy of the 
localized particles and the last term accounts for the in- 
teractions due to the planar displacements ui = {u X i, u V i) 
of the electrons around their equilibrium positions: 
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Expanding the last term for small displacements re- 
sults in a series expansion for the energy Eq. Q in 



powers of l/rl^ 2 The leading term, proportional to 
l/r s , corresponds to the Madelung energy E M of Eq.(J2)l. 
In free space, it attains its minimum value Ebcc — 
— 0.89593/r s (in atomic units) for a Body Centered Cubic 
(BCC) Wigner crystal^ The second term in the expan- 
sion, proportional to \/r Z J 2 , is the zero point energy of 
the particle fluctuations in the harmonic approximation, 
which also depends on the selected crystal structure. It 
is negligible for r s — > oo, and remains smaller than the 
Madelung term by typically an order of magnitude at 
r s ~ 100. Nonetheless, it can play an important role in 
determining the relative stability of the different crystal 
structures, especially when approaching the melting den- 
sity. Higher orders in the energy expansionist include 
anharmonic (l/r s p with p > 2) and exchange terms of 
the form e~ c ^, which we shall neglect in the following 
discussion. 

Up to quadratic order in the displacements, our model 
Hamiltonian reads: 
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where Zy is a 2 x 2 matrix characterizing the dipole-dipole 
interactions, given by (a, (3 — x,y): 
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with Rij = Ri—Rj. The most general elementary Bravais 
lattice compatible with a given layered structure is iden- 
tified by a couple of basis vectors describing the ordering 
within the planes, A\ — (a±,0, 0), A2 = (o^, a2y, 0), 
and a third vector A3 = (03^ , a^ y , d) which sets the rel- 
ative shift (a3 X ,a3y) between two equivalent 2D-lattices 
on neighboring planes. Other structures, with more than 
one particle per unit cell, are possible in principle, but 
will not be considered here. 

Due to the additional lengthscale d introduced by the 
layered constraint, the crystal energy is no longer a func- 
tion of r s alone. Its dependence on the lattice geometry is 
best expressed by introducing a dimensionless parameter 
7, which measures of the relative importance of interlayer 
and intralayer interactions. It is defined as the ratio be- 
tween the mean interparticle distance in the planes and 
the interlayer separation, namely 7 = ypKr s ,iD I 'd. Here 
r s ,2D defines the 2D density parameter in the individual 
layers, related to the bulk r s by r 2 s2D = 4rj?/3d. The 
first two terms of the low-density expansion, correspond- 
ing respectively to the Madelung energy and the zero- 
point fluctuation energy in the quadratic model 10} can 
be written in compact form as: 
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It should be noted that an effective mass m* ^ m and 
a dielectric constant k / 1 can be straightforwardly in- 
cluded in the model through a redefinition of the Bohr 
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radius as 



aBn{m/m*), unit energy me A /K 2 



m*e /K 2 h 2 , and density parameter r s — ► r s {m* /m)/ n. 
Hereafter, energies and lengths will therefore be ex- 
pressed in terms of these effective units, characterizing 
the host medium. A much more complex situation arises 
in systems with a frequency-dependent dielectric screen- 
ing, leading to the formation of polarons, for which the 
reader is referred to Refsii£*2£. 



B. Minimization of the Madelung energy 

Following the hierarchy of the series expansion intro- 
duced above, we start by searching for the layered config- 
uration which minimizes the electrostatic repulsion be- 
tween the particles, which is appropriate in the limit 
of large r s . The calculation is performed by standard 
Ewald summation techniques, which split the slowly con- 
vergent series in Eq. J3J) into two exponentially converg- 
ing sumsS Given the interlayer separation d and the 
bulk density n (or, alternatively, given the pair of di- 
mensionless parameters 7 and r s ) we are left with 4 free 
minimization parameters: 2 for the inplane structure, 2 
for the interlayer ordering. 

The result of the minimization for the Madelung co- 
efficient A in the range < 7 < 6 is illustrated in Fig. 
^ Two distinct regimes can be identified. In the limit of 
large separations (7 < 1), the coupling between the lay- 
ers is weak, and the resulting planar pattern is hexagonal, 
with a staggered interlayer ordering, i.e. the particles on 
the neighboring layers falling on top of the centers of 
the triangles. The sharp rise of the Madelung constant 
in this regime is due to the fact that the compensat- 
ing background is distributed homogeneously in three- 
dimensional space, which penalizes strongly anisotropic 
charge distributions^ 

Upon reducing the interlayer separation so that 7 > 1, 
the increasing interlayer interactions make the hexago- 
nal pattern energetically unfavorable. Above 7 = 1.15, 
a more isotropic ordering of the charges is stabilized, 
which presents a centered rectangular (CR) structure in 
the planes. Further increasing 7 leads to a sequence of 
structures whose planar patterns are respectively squared 
(S, in the interval 1.32 < 7 < 2.13), rectangular (R, 
2.13 < 7 < 2.84), centered rectangular (CR, 2.86 < 
7 < 4.31), a generic rhombic, or oblique phase (Rh, 
4.31 < 7 < 4.45), then hexagonal again, and so on. Such 
phases are all connected by continuous structural tran- 
sitions, with the exception of the hexagonal structure, 
which is attained through a discontinuous change of the 
crystal parameters. Note that in the very narrow interval 
2.84 < 7 < 2.86, a generic structure with rhombic planar 
symmetry is stabilized, which allows to evolve continu- 
ously from the rectangular to the centered rectangular 
patterns (not shown) . The sequence of structural transi- 
tions goes on at larger values of 7. 

The interlayer ordering is shown at the bottom of Fig^ 
It is staggered for the first three patterns (H), (CR) and 
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FIG. 1: Madelung coefficient A in atomic units for differ- 
ent crystal structures constrained to a layered environment, 
as a function of the anisotropy ratio 7. The different curves 
correspond to different planar configurations: hexagonal (H), 
square (S), centered rectangular (CR), rectangular (R) and 
rhombic (Rh). The interlayer orderings in the simplest cases 
at low 7 are sketched below the curves (for the R and CR 
structures, the stacking varies as indicated by the double ar- 
rows). The resulting three-dimensional Wigner crystal re- 
duces to a perfect BCC at the points marked by filled dots, 
whose energy is indicated by the horizontal arrow. 



(S) for 7 < 2, as expected for large interlayer separa- 
tions, where the relative ordering is fully determined by 
the coupling between two adjacent planes, and indeed 
coincides with what is found in bilayer systems 1 ^ (see 
Section III El below h At larger values of 7, the interac- 
tions beyond the nearest planes become relevant, which 
makes the simple staggered ordering unfavorable. For in- 
stance, a staggered/non-staggered transition takes place 
within the rectangular phase at 7 = 2.46, corresponding 
to a relative sliding of the planar structures on adjacent 
planes in the direction of the long bonds (indicated by 
the double arrow in Fig. 1). 

Remarkably, each of the phases identified above con- 
tains a special point 7* where the ideal BCC structure — 
which has the lowest possible Madelung energy in three 
dimensions — is itself compatible with the layered con- 
straint. The different planar configurations identified 
above then correspond to the different ways of cutting 
a BCC by an array of equally spaced layers. Such points 
are easily calculated by setting the distance d = 2n/\K\, 
with K any reciprocal lattice vector, and correspond to 
7* = 2 1 / 4 , 2, 2 1 /4 3 3/4 ;2 i/4 5 3/4 j 2 33/4^ etc. .. Similarly, 

the higher relative minima visible in Fig. ^ correspond 
to different orientations of the same three-dimensional 
Face Centered Cubic (FCC) ordering. 

Away from such special points, the overall charge dis- 
tribution remains very isotropic in all the region 7 > 1, 
as testified by the extremely small deviations of the 
Madelung energy from the ideal case, AEm < 10~ 4 /?v 
Such small energy variations, however, refer to the op- 
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timal structures obtained at different vaiues of 7, which 
does not mean that the electrostatic repulsion between 
the carriers is irrelevant in the determination of the 
charge ordering patterns in real systems: in a given 
compound, where both the interlayer distance and the 
density arc fixed, one should rather compare the ener- 
gies of two competing phases at fixed 7. For exam- 
ple, enforcing a hexagonal symmetry at 7 = 2, where 
the optimal structure is squared, would cost an energy 
AE M ~ 0.015/r s ~ 2QQK at r s = 20, which is compa- 
rable with the typical charge ordering energy scales in 
solids. Yet, since the Madelung energy is determined by 
the interactions with a large number of (distant) neigh- 
bors, the structures found here are expected to be rel- 
atively soft against local deformations. The situation is 
different regarding global symmetry changes, as can re- 
sult from the inclusion of a periodic potential of compet- 
ing symmetry, which could strongly modify the sequence 
and order of the structural transitions, possibly favoring 
the appearance of alternative phasesiS^S 



C. Zero point fluctuation energy 

The next term in the series expansion of the ground 
state energy Eq. © corresponds to the quantum zero 
point fluctuations of the particles around their equilib- 
rium positions, in the harmonic approximation. It is 
negligible at large r s (low density) , but it becomes quan- 
titatively important at lower r s , where it can slightly 
modify the sequence of phases identified in the preceding 
Section. Upon further reducing r s , this term eventually 
drives the quantum melting of the crystal, that will be 
analyzed in the next Section. 

The calculation of the fluctuation term proceeds as fol- 
lows. The harmonic model Eq.Q) is diagonalized by in- 
troducing the normal modes q s ^ 
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FIG. 2: a) Zero point vibrational term B for the different 
structures identified in Fig. 1, within their ranges of me- 
chanical stability. The sequence of structures with the lowest 
vibrational energy is indicated at the bottom. The arrow 
indicates the value (2/3)B (3D) = 0.887, where B ( ' iD) is the 
vibrational energy of a BCC crystal in vacuum; b) Inverse 
moment M_i of the DOS, which is proportional to the mean 
electronic fluctuation (u 2 ) Eg. 1121 for the same structures. 
The arrow indicates the value for a BCC in free space. 
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where s s ? are the two-dimensional polarization vectors 
(the electrons oscillate within the planes) and the vec- 
tor k runs through the Brillouin zone of the three- 
dimensional reciprocal lattice. This yields two branches 
s = 1, 2 of collective modes with eigenfrequencies tu s so 
that the vibrational energy per particle can be expressed 
as: 
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It is useful to introduce the normalized density of states 
(DOS) of the collective modes, that we write here in gen- 
eral as: 
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(D is the number of branches, corresponding to the di- 
mensionality of the electron motion) as well its dimen- 
sionless moments: 



M n = Jduj p(u>) (j^- 



(10) 



with ujp = 3e 2 /(mri?) the usual 3D plasma frequency. 
With these definitions, the vibrational energy in Eq.|0 
is seen to be directly proportional to the first moment of 
the DOS, with 
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The usual 3D case in vacuum is recovered by restoring 
the out-of-plane oscillations in Eq. (4), and by setting 
D = 3 in Eq. (|11|) . For example, for the BCC structure 



we find M 



(3D) 



0.511, which yields the well known value 



= 1.33M 
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The analysis of the frequency spectrum shows that 
each given structure has a limited interval of mechani- 
cal stability: for certain geometries, the dynamical ma- 
trix acquires negative eigenvalues around some critical 
wavevector k c , corresponding to purely imaginary col- 
lective frequencies which drive the crystal unstable (this 
phenomenon also exists in free space, where FCC and 
the simple cubic structure are known to be intrinsically 
unstable). For example, in the interval of 7 under study, 
a structure with hexagonal symmetry is only stable for 
7 < 1.32, 3.5 < 7 < 4.95 and 5.05 < 7 < 5.8. 

We have calculated the fluctuation term B(j) for the 
different symmetric structures (H, R, S, CR) identified 
in the previous section, within their respective intervals 
of mechanical stability, as well as for the rhombic phase 
at 2.84 < 7 < 2.86 and 7 > 4.31, which is shown in 
Fig. 0a. As for the Madelung energy, two essentially 
different regimes can be identified. For 7 < 1, the elec- 
tron motion is mostly determined by the Coulomb in- 
teractions within the layers (interlayer forces are neg- 
ligible) and the collective modes of the pure 2D case 
are recovered. If normalized by an appropriate "two- 
dimensional plasma frequency" w\ D — e 2 /mrj? 2D , the 
first moment in the hexagonal phase tends to the con- 
stant value Mi y2 D = 0.814ml Going back to the present 
three-dimensional units, however, where the moments are 
normalized as in Ea. (|10f> . the fluctuation term diverges 
at large separations as ^(7) ~ 7r 1 / 4 3 1 / 2 Mi^d/^1 1 ^ 2 ■ In 
the regime 7 > 1, on the other hand, the fluctuation term 
flattens around a value which roughly corresponds to 2/3 
of the fluctuation in free space, indicated by the arrow in 
Fig. 12 a. This follows from the fact that only the oscilla- 
tions along 2 of the 3 space directions are allowed, as we 
can see explicitely from Eq. 1)110. 

The structural phase diagram resulting from the analy- 
sis of the total energy (jSJ, including the vibrational term 
(|11[) . and taking into account the ranges of mechanical 
stability of the different phases, is reported in Fig. [3J The 
first observation is that, apart from the disappearance of 
the CR phase from certain intervals, which is penalized 
by its higher vibrational energy than the H phase, the 
locus of the structural transitions does not change much 
with r s . The sequence of phases identified in Fig. 1, 
based on the analysis of the Madelung energy, is recov- 
ered at extremely large values of r s . On the other hand, 
the vibrational term affects the structural transitions al- 
ready at r s < 1000. This is due to the fact that, even 
though the electrostatic term A/r s is still larger than the 

zero-point fluctuation energy B/r^ 2 , the latter under- 
goes much larger relative variations among the different 
phases. Below r s ~ 100, the phase diagram is entirely de- 
termined by the minimization of the vibrational energy 
(see Fig. 0a). As was stated above, however, the overall 
shape of the phase diagram does not depend much on r s , 
the transitions being essentially determined by the pa- 
rameter 7. Let us also remark that the vibrational term 
is much less influenced than the Madelung term by the 
specific interlayer arrangements, whose effect (if any) is 
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FIG. 3: Structural phase diagram of the Wigner crystal in 
a layered environment, based on the total energy as a 
function of the anisotropy ratio 7 and the bulk density pa- 
rameter r 3 . The labels are the same as in previous figures. 
The solid (dashed) lines are for structural transitions where 
the crystal parameters evolve discontinuously (continuously) . 
The bold lines indicate mechanical instabilities, accompanied 
by a discontinuity of the crystal energy. The melting line is 
determined by solving Eq. 1131 . For the hatched region, see 
text. 

to slightly modify the range of mechanical stability of 
each phase. 

Another fundamental property of the system, which 
gives valuable informations on the collective vibrations 
of the particles, is the mean electronic fluctuation (u 2 ). 
In the harmonic approximation, this quantity is propor- 
tional to the inverse moment of the DOS of the collective 
modes, defined in Eq. i|10|): 

(u 2 ) = — S T- = T^-rl' 2 12 

where again we keep track of the explicit dependence 
on the dimensionality D. As can be seen in Fig. |2Jb, 
it increases as each phase approaches the boundaries of 
its stability range. This is because the mechanical in- 
stabilities are approached via a softening of a branch of 
collective modes, causing an increase of the DOS at low 
frequency and, through Eq. I|10|) . of the inverse moment 
M_i. A local increase also occurs at the points where the 
staggered interlayer ordering is lost (see e.g. the maxi- 
mum at 7 = 2.46 within the R phase in Fig. Elb). 

From analogous arguments, it follows from Eqs. i|10|) 
and (|11|) that the vibrational energy generally attains its 
minimum value close to mechanical instabilities. Within 
the present approximate framework, this can cause the 
total energy to jump discontinuously at the instability 
point when the next stable phase is attained, which cor- 
responds to the bold lines is Fig. 3. For example, the 
hexagonal lattice becomes unstable at 7 > 1.33 and, for 
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r s < 100, the transition to the square phase is accompa- 
nied by a small jump in energy. Such discontinuities can 
in principle be avoided by allowing for Bravais lattices 
with more than one electron per unit cell (the resulting 
internal structure could then be assimilated to some lo- 
cal tendency to electron pairin g??'? 4 ). Note also that it 
is precisely close to mechanical instabilities, where (u 2 ) 
is largest, that the neglected anharmonic corrections are 
expected to be most important. Their consequences on 
the structural phase diagram presented here deserve fur- 
ther theoretical study. 



D. Melting of the crystallized state 

In this section, we analyze the melting of the crystal- 
lized state by making use of the Lindemann criterion, ac- 
cording to which a transition to a liquid phase takes place 
when the spread (u 2 ) attains some given fraction S of the 
nearest-neighbor distance a n . n .. We take 8 = 0.28 from 
Ref^, which is appropriate for the quantum melting of 
b oth 2 D and 3D Wigner crystals. Solving the equation 
\J (u 2 ) /a n , n , = 5 in terms of the density parameter r S) 2D 
in the planes leads to: 



' s,2D 



M_i( 7 ) 1/2 
28 2 C 2 { 1 ) 



(13) 



where C — a n , n ,/r s ,2D is an aspect ratio relating the 
nearest-neighbor distance to the density parameter in the 
planes, and the implicit condition 7 = \pKr c s 2D jd holds. 
Note that for structures with rectangular symmetry, the 
Lindemann criterion must be modified to account for the 
existence of two nonequivalent near-neighbor distances. 
Here we use a simple generalization which consists in 
replacing a n . n . with the average of the two shortest near- 
neighbor distances, and which reduces to the ordinary 
criterion for the square and hexagonal structures. A 
check of the validity of such generalized Lindemann cri- 
terion will be given in Section II E, by direct comparison 
with independent theoretical results on bilayer systems. 

The melting curve deduced from Eq. Ijl3(l for the differ- 
ent structures considered here is illustrated in Figs [21 and 
The most important result is that the crystal melt- 
ing can be pushed to higher densities by reducing the 
interlayer spacing, which can already be inferred by ne- 
glecting the weak 7-dependence of the coefficients C and 
M_i of Eq. (|13l) in the region 7 > 1. The main reason 
to this is that for 7 > 1 the electron spread is essentially 
governed by three-dimensional Coulomb interactions, as 
we can see from the explicit dependence of Eq. (|12|) on 
the bulk r s , while the electron motion is two-dimensional, 
so that the appropriate nearest-neighbor distance for the 
Lindemann ratio is proportional to the planar density 
parameter r s ^D — (2/V3d)r^ 2 m& 

In addition, for each given spacing d, the geometrical 
confinement leads to a further stabilization of the crystal 
through a reduction of the spread (u 2 ) itself. This effect 



is directly reflected in Fig. [5]b in a reduced value of M_i 
as compared to the corresponding value in free space, and 
should not be confused with the trivial dimensional factor 
D, that was taken out explicitely from Eq. (|12ll . It is due 
to the fact that, as soon as the cubic symmetry is lost, the 
restoring forces induced by the dipole-dipole interactions 
Eq. (5) are not equivalent in the three space directions, 
so that the electron fluctuation becomes anisotropic (the 
observed shrinking of the planar spread would occur at 
the expense of increasing the out-of-plane fluctuations, 
which are anyhow suppressed in the model). To give an 
example, taking an average value M_i ~ 2.4 and the as- 
pect ratio C = \pH for the square planar ordering yields a 
critical value r c s 2D ~ 4.9%/d. Comparable results (within 
few percent) are found for the other structures. 
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FIG. 4: Phase diagram of the Wigner crystal in a layered envi- 
ronment, as a function of the interlayer distance d (in units of 
the effective Bohr radius a* B ). The discontinuity close to the 
three-phase critical point (S, H, liquid) is due to the different 
aspect ratios C of the two competing structures. The hatched 
area is a region possibly characterized by an anisotropic liquid 
behavior (see text). The shaded area corresponds to 7 > 6 
and has not been studied. 



In the opposite limit of large separations (7 <C 1), 
where interlayer forces become negligible, we recover the 
usual critical value r c s 2D ~ 40 for the 2D hexagonal 
Wigner crystal. Note that the actual critical value at 
finite 7 always lies below this asymptotic estimate, con- 
firming that the inclusion of interlayer interactions causes 
a stabilization of the crystal phase compared to the pure 
two-dimensional case, as was argued in the introduction. 

A few comments on the limits of validity of the present 
model are in order. First, the enhancement of Wigner 
crystallization predicted by Eq. I|13|l cannot extend indef- 
initely: the melting line should eventually saturate at low 
separations when isotropic electron motion and three- 
dimensional screening are restored by interlayer tunnel- 
ing processes^ On the other hand, as was stated in the 
introduction, replacing the host lattice of ions by an ef- 
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FIG. 5: Lindemann ratio y/ (u 2 )/a n , n , as a function of 7 
(r s ,2D, upper abscissa) at a given interlayer spacing d = 20. 
The continuous line is the average Lindemann ratio, the 
dashed line represents the Lindemann ratio along the direc- 
tion of the closest near-neighbor (see text). The horizontal 
line sets the critical value for melting. Upon increasing the 
electron density, the transition from the crystal to the liquid 
could occur through an intermediate anisotropic liquid phase 
(hatched region, see also Fig. 4). 
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FIG. 6: Phase diagram for the symmetric electron bilayer, 
in terms of the two dimensional density parameter r s .2o, as 
a function of the interlayer spacing d. The lengths are scaled 
to atomic units. H,S,R denote respectively the Hexagonal, 
Square , Rectangular phase. Note that the QMC simulation 
of Ref.— ^ was restricted to study only two phases (H,S), while 
an additional rhombic phase (Rh) could be stabilized in the 
DFT calculations of Refill. 



fective jellium is allowed provided that the spread of the 
electron wavefunction is larger than the ion-ion distance 
<iq. From eq. l|12|) . the condition \J (u 2 ) > ao gives r s > 4 
{ r s.2D ^ 6) for a typical value of ao = 3 A, assuming k = 1 
and m* — m. Below this value, the discrete nature of the 
host lattice should be included, which can further stabi- 
lize the crystallized state, as pointed out in Refs>i 7 ^i£*2&. 

Before concluding this section, let us remark that, for 
anisotropic planar orderings such as the rectangular and 
the centered rectangular structure, two independent Lin- 
demann ratios could in principle be defined (one for each 
nonequivalent near-neighbor direction) rather than the 
single average criterion used so far. It would then appear 
that the melting along the short bonds is much easier 
than along the long bonds, due to the closer overlap be- 
tween the electron wavefunctions. This phenomenon is 
illustrated in Fig. [S] and could imply a tendency towards 
an anisotropic (or "striped") liquid phase, which is gener- 
ally not ruled out by the isotropic nature of the Coulomb 
repulsion (see also the hatched regions in Figs. 3 and 
4) ,37 rp^g resu ^ s reported in Fig. 5 also indicate a pos- 
sible reentrant behavior, although no conclusive answer 
can be given at this level of approximation. 



E. Symmetric electron bilayer 

In this section we analyze a system composed of two 
coupled electronic layers, in order to check the validity of 
our approach by direct comparison with available Density 
Functional TheorjiSi and Quantum Monte Carlo based 
calculations^^. In the early work on classical bilayersji^ 



the analysis of the Madelung energy showed that several 
structural phase transitions occur as the distance d be- 
tween the two planes is varied while keeping the electron 
density fixed. At short distances the planar ordering is 
rectangular, and collapses to the usual hexagonal phase 
in the formal limit d — > 0. This phase evolves continu- 
ously into a staggered square structure, when d is of the 
order of the interparticle spacing, which is clearly remi- 
niscent of the BCC structure observed in 3D space (cf. 
the discussion in Section II B). Upon further increasing 
d, the lattice progressively deforms into a rhombic phase, 
to attain the hexagonal staggered phase expected in the 
limit of independent layers. Including the zero-point en- 
ergy of the collective excitations as in Eq. © raises the 
energy of the rhombic phase, which therefore disappears 
from the phase diagram at sufficiently high density, leav- 
ing the other transitions essentially unchanged. 

We have analyzed the quantum melting of the different 
Wigner crystal structures realized in such bilayer system 
by making use of the Lindcman criterion discussed in 
the preceding Section. We see from Fig. 6 that both 
the sequence of phases and the critical melting densi- 
ties obtained within the present quadratic approxima- 
tion are in satisfactory agreement with the more sophis- 
ticated numerical results of RefsiSiiSSiSi (the melting den- 
sity is slightly underestimated as compared with QMC, 
but quite similar to the DFT result). It is interesting to 
see that the same trends observed in the preceding Sec- 
tion for the layered solids are already present in the single 
bilayer. In particular, reducing the interlayer separation 
leads to a sensible stabilization of the crystal compared 
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to the isolated layers. This is clear in Fig. where 
the the melting line always lies below the critical value 
r c s 2D ~ 40 of a purely 2D Wigner crystal. Note also 
that, contrary to Refi22, we find that the enhancement 
of Wigner crystallization is slightly more pronounced in 
an infinite array of layers than in a single electron bilayer. 



III. WIGNER CRYSTALLIZATION IN QUASI 
ONE-DIMENSIONAL SOLIDS 



We now extend our analysis to the case of quasi one- 
dimensional solids, which we model as periodic arrays 
of conducting wires. Following the general arguments 
presented in the previous Section, the enhancement of 
Wigner crystallization in this case should be even more 
pronounced than in the two-dimensional case, because 
of the suppression of electronic motion in two transverse 
directions rather than one. The effect is even more dra- 
matic if we consider that a quantum crystal with gen- 
uine long-range order cannot be realized in a pure one- 
dimensional system^ while it is stabilized if we account 
for the long-range Coulomb interactions between carriers 
on different wires 44 

We shall consider here a square array of wires for il- 
lustrative purposes, although the specific arrangements 
occurring in real solids (rectangular, rhombic) can be 
treated case by case. Assuming a simple ordering of pe- 
riod a within the wires and an interwire distance d, the 
most general elementary three-dimensional Bravais lat- 
tice compatible with the given geometrical constraint is 
described by the following basis vectors: A\ = (0,0, a), 
A 2 = (d, 0, b), A 3 = (0, d, c). The volume of the 3D uni- 
tary cell is V c = ad 2 = 47rr 3 /3, the anisotropy ratio is 
now defined as 7 = a/d and the ID density parameter is 
r s .iD = a/2. As in the layered case, we take a compen- 
sating positive charge distributed uniformly in the bulk. 
The analysis presented in the preceding Section can be 
repeated here following the same steps: i) calculation 
of the structure with the lowest Madelung energy upon 
varying the anisotropy ratio; ii) calculation of the corre- 
sponding vibrational energies; iii) determination of the 
melting curve via the Lindemann criterion. The gener- 
alization is straightforward, and we only report here the 
main results. 

The stuctural phase diagram (Fig. |JJ is clearly less rich 
than in the layered case, because once the density and the 
interwire distance d are fixed, only the relative ordering 
between the electronic crystals on neighboring wires re- 
mains to be determined, corresponding to the pair of pa- 
rameters b and c. In the limit 7 — > 0, the interwire inter- 
actions vanish and the limit of isolated wires is recovered: 
the Madelung constant A diverges due to the isotropic 
distribution of the jellium, as explained previously (cf. 
footnote 43). In this limit the interwire ordering is stag- 
gered, with b/a — c/a = 1/2, corresponding to a body 
centered tetragonal (BCT) lattice in three-dimensional 
space. The BCT structure, everywhere compatible with 
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FIG. 7: Phase diagram for a three-dimensional Wigner crystal 
embedded in a square array of 1-dimensional wires, of side d. 
Lengths are scaled to effective atomic units. For the definition 
of the phases I and II, see text. The shaded region corresponds 
to 7 > 10 (not studied). 



a square array of wires, has the lowest Madelung energy 
in the whole range < 7 < 2.83, with the two special 
values 7* = v2 and 7* = 2 corresponding respectively 
to a BCC and a FCC. For 7 > 2.83 the minimum con- 
figuration becomes less symmetric, with c/a ^ 1/2 but 
the ratio b/a still locked to the value 1/2 up to 7 = 6.99. 
This phase is denoted (I) in Fig[7| Beyond 7 = 6.99, a 
second structural transition occurs leading to a generic 
phase (II) with both b/a ^1/2 and c/a ^ 1/2. Other 
transitions can take place at larger values of 7, within the 
generic phase II. The sequence of phases does not change 
upon inclusion of the vibrational term. 

By applying the Lindemann rule we obtain a paramet- 
ric formula for the melting curve analogous to Eq. 1131) : 



1 



' s,lD 



(128tt) 1 /3 



M_!( 7 ) 

6 2 



2/3 



d 2 ' 3 



(14) 



with the implicit condition 7 = r c s 1D /2d. The conse- 
quences of geometrical confinement evidenced in the lay- 
ered case are recovered here. The electron spread along 
the wires is again governed by the three-dimensional 
plasma frequency [cf. Eq. I|12|) ]. due to the isotropic 
nature of the Coulomb interactions, while the nearest- 
neighbor distance here scales with r s ,iD = / 3)r 3 / d 2 . 
Further stabilization of the crystallized state is achieved 
through a reduction of the electron spread along the 
wires, revealed by an inverse moment M_i which is 
typically 50% lower than the value in vacuum. Its 7- 
dependence for 7 > 1 is quite flat (not shown), except 
in the vicinity of the transition at 7 = 2.85, where it 
raises due to the mode softening discussed in Section II 
C. Replacing the average value M_i ~ 2 into Eq. i|14|) 



yields 



s.lD 



1.2c? 2 / 3 , corresponding to an even stronger 
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enhancement of Wigner crystallization than in the lay- 
ered case (see Table QJ. In the opposite anisotropic limit 
7 « 1, M_i diverges as in the case of an isolated wire 
(cf. footnote 43), so that the Wigner crystal is never 
stabilized (r c s 1D — > oo). 
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TABLE I: Definition of the anisotropy ratio 7, approximate 
melting lines obtained for quasi two-dimensional and quasi 
one- dimensional systems, and specific values obtained at two 
different interlayer (interwire) distances d, expressed in units 
of the effective Bohr radius a* B (right columns). 



IV. CONCLUSIONS 

We have investigated the Wigner crystallization of 
electrons in quasi low-dimensional compounds, where 
the carrier motion is effectively low-dimensional, while 
the Coulomb interactions are assumed long-ranged and 
isotropic. The system properties are found to depend 
crucially on the ratio 7 of the mean interparticle spacing 
within the conducting units (layers or chains) to the sep- 
aration d between units. While the behavior expected for 



isolated units is recovered at large separations (7 <C 1), 
an overall isotropic ordering of the charges is achieved 
for 7 > 1, when the interactions between different units 
become important. In this case, three-dimensional struc- 
tures as close as possible to the ideal case of a BCC 
are formed, leading to a cascade of structural transitions 
which can be tuned by varying the particle density, or 
the distance d itself. In addition to this rich phase di- 
agram, the presence of isotropic Coulomb interactions 
in such anisotropic compounds results in a strong sta- 
bilization of the charge ordered phases, possibly up to 
densities of practical interest, where the characteristic 
energy scales of the Wigner crystal can become compa- 
rable with other relevant scales in the solid. Although it 
is clear that the interplay with several other factors such 
as the periodic lattice potential^*iSiiLi2i2S^ chemical 
impurities^ polaronsi^Si or magnetic interactions^^ 
should be considered for an accurate description of real 
materials, the long-range Coulomb interactions appear 
in light of the present study as a key ingredient to un- 
derstand the charge ordering phenomena in quasi low- 
dimensional systems. 
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The kinetic energy evaluated in the noninteracting limit is 
strongly reduced in lower dimensions, where there are less 



high momentum states available (T = 1.11/ r 2 , 0-5 /r 2 t2 D 
and 0.11/rf 1D in 3D, 2D and ID respectively) while the 
Madelung energy in the opposite classical limit is less de- 
pendent on dimensionality [Ezd = 0.896/r s for the body 
centered cubic lattice, E^d = l.ll/r s ,2r) for the hexagonal 
lattice]. Incidentally, the crystallization transition in both 
3D [r s — 100) and 2D (r s ,2D — 40) takes place at the same 
value of the ratio E/T ~ 81. 

The usual value for the Madelung energy of the hexagonal 
lattice in two dimensions would be recovered at large in- 
terlayer separations by considering a layered compensating 
charge. This would not alter the sequence of crystal struc- 
tures, since it would add a constant term to the energy 
depending only on the anisotropy ratio 7. 
This can be understood by observing that the logarith- 
mic divergence of the mean spread (u 2 ) ~ J dwuj" 1 arising 
in ID is washed out in the present case where the wires 
are embedded in a three-dimensional environment and mo- 
mentum integrations run over the 3D Brillouin zone. 



